Coexistence of exponentially many chaotic spin-glass attractors 
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A chaotic network of size TV with delayed interactions which resembles a pseudo- inverse associative 
memory neural network is investigated. For a load a — P/N < 1, where P stands for the number 
of stored patterns, the chaotic network functions as an associative memory of 2P attractors with 
macroscopic basin of attractions which decrease with a. At finite a, a chaotic spin glass phase 
exists, where the number of distinct chaotic attractors scales exponentially with A'^. Each attractor 
is characterized by a coexistence of chaotic behavior and freezing of each one of the A'' chaotic units 
or freezing with respect to the P patterns. Results are supported by large scale simulations of 
networks composed of Bernoulli map units and Mackey-Glass time delay differential equations. 



Over two decades ago a link between statistical physics 
and neural network theory was established, where as- 
sociative memory models were mapped onto magnetic 
Hamiltonian systems with random and frustrated inter- 
actions [ll, m . The capacity, i.e. the number of stored 
patterns, was calculated using techniques borrowed from 
random frustrated magnetic systems and in particular 
advanced method to solve infinite-range spin-glass mod- 
els 0. In parallel to the progress in interdisciplinary 
research on random magnetic systems, the theory of non- 
linear dynamics developed and recently there have been 
important steps forward in the understanding of condi- 
tions to achieve chaotic networks synchronization 0, ■ 

It is reasonable to assume that the two phenomena, 
spin-glasses and chaos, represent two conflicting trends. 
The spin-glass (SG) phenomenon is characterized by the 
existence of exponential number of meta-stable states, 
in which the local degrees of freedom (spins) are frozen 
0, S 01- C)n the other hand, chaotic dynamics is char- 
acterized such that two nearby trajectories diverge from 
each other as a consequence of the positive Lyapunov ex- 
ponent [1]. Nevertheless, a network of coupled chaotic 
units can synchronize 01 @ to an identical chaotic tra- 
jectory. In this paper we show that chaotic networks can 
have a SG phase, characterized by an exponential number 
of chaotic attractors with the system size, N. Each one 
of the exponentially many attractors is characterized by 
a coexistence of chaotic behavior and substantial freezing 
of each one of the N chaotic units. We start in section |T] 
with the definition of the model for a network of chaotic 
maps. The associative memory features of the model are 
analyzed at section [TT] and the spin glass phase is exam- 
ined at section imi The generalization of the chaotic map 
networks to time delay differential equations is discussed 
in section [TVl and concluding remarks are given at section 
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I. MODEL 

One of the first examined neural network models 
within the statistical mechanics framework is an asso- 
ciative recall of memory without error known also as 



pseudo-inverse model 0, [T^] ■ The model consists of a 
fully connected network governed by the Hamiltonian 



H 
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N 



(1) 



where S'i = ±1 is an Ising spin and the sum is over all 
pairs of spins. The coupling matrix Wij is defined as 



I P _i , 1 ^ 



(2) 

where N is the size of the network. The stored patterns 
= 1 . . . P, i = 1 ■ ■ ■ N) are P vectors of length N 
and = ±1 with equal probability. All the patterns 
are minima of the Hamiltonian and the pseudo-inverse 
network operates as an associative memory up to a 
P/N = l [lOf, which is the upper bound for the capacity 
of such networks with symmetric interactions. In the case 
where the diagonal interactions, Wu, are introduced the 



capacity was found to be a = 0.5 [9|. The matrix Wij 
has P eigenvalues equal to 1 where the patterns are the 
corresponding eigenvectors. 
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and all other N — P eigenvalues are identically zero. 

The counterpart chaotic map version of the pseudo- 
inverse model is given for unit i, for instance, by 

N 

Xl = (1 - e)HXl-') + eY^W,,g{X]--) (4) 

where Wij is given by ^ including the diagonal ele- 
ments, e is the strength of the coupling, 1 — e is the 
weight of the internal dynamics and and Q{-) stand 
for (non-linear) response functions. For simplicity of the 
discussion below we mainly present results for the sim- 
plest chaotic map, the generalized Bernoulli map [llj 



^{x) ^ Q{x) ~ sign{x) (^{a\x\) mod l) 
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FIG. 1: Raster diagrams of the jumps during time evolution 
of a network (A'' — 40, P = 20) of chaotic BernoulU maps, 
eqs. (|4l5p . (a) An initial condition has a unity overlap with 
one pattern only, eq. ((Tjl. Jumps occur simultaneously for 
all units and the overlap with the pattern remains 1. (b) 
Random initial condition. Units jump asynchronously and at 
a different rate and the maximal overlap with the patterns is 
less than 3/y/N. Note that in case (b) there are more jumps as 
time evolves since the units get frozen, mostly, in the interval 
Xi e [0.5,1] (section 11111. 




FIG. 2: (color online). Estimation of the basin of attrac- 
tion. The average minimal initial overlap with one of the 
patterns which converges with probability 0.9 to an overlap 
> 0.9 with the same pattern, iV = 20(0), 50(o), lOO(x). Sim- 
ulation parameters are: r = 40, a — 1.1, e = 0.9, each mea- 
sured point was averaged over at least 10^ initial overlaps, 
the system evolved for lO^r steps before overlaps were aver- 
aged over a window of r steps. Inset: Measured Lyapunov 
exponent. A, for both Mattis (an overlap close to unity with 
one of the pattern) and SG chaotic attractors, normalized to 
Aoo = r-Mn|ae/(l - a(l - e))| for r = 400, iV = 50 and each 
point averaged over 50 attractors. 



where J^(a;) G [^1,1] and the map is known to be chaotic 
for a > 1. For all simulations exemplified below the se- 
lected parameters are a = 1.1, t = 40 and a = 0.5, 
unless otherwise mentioned. 

The state vector X* of the N miits at time t can be 
written as 

p 

where S]_ represents the orthogonal part of X* to the sub- 
space spanned by the P(< N) patterns. Since Wij has 
zero eigenvalues in directions perpendicular to the sub- 
space spanned by the patterns, decays exponentially 
to zero for the discussed region a(l — e) < 1. 

II. MATTIS STATES AND ASSOCIATIVE 
MEMORY 

An initial condition equals to one of the patterns is a fix 
point of the Hamiltonian system eqs. (1-2). Nevertheless, 
for the chaotic dynamics, eqs. (j4l5p . an initial condition 
which has a unity overlap with pattern /i 

X* = b^^^ {l^t^T, hi > 0) , (7) 

does not lead to a fixed point. Using the fact that 

HKO = C-^(^m) and Ef=i W,,C^ = equation © 
is reduced to 

bl^il~e)^ibl~') + enb'r) (8) 



which is known to be chaotic for a > 1 and describes 
merely the behavior of the coefficient b^. For binary pat- 
terns, ilj jumps of the Bernoulli map (i.e. \Xf \ > 1/a) 
occur for all units simultaneously as exemplified by Fig. 
[T^. The network has has solely a unity projection on 
pattern /x, and the N dimensional dynamics remains in 
the direction of the pattern. 

For a random initial condition of the system 
composed of projections to all patterns, X* = 
b^i^ (1 '^t), jumps do not occur simultane- 
ously for all units. Fig. [TJd. The dynamics of the coeffi- 
cients {6^} are coupled and cannot be described indepen- 
dently as in ([8]). A nontrivial effective coupling among 
the patterns is generated and time dependent overlaps 
with the P patterns is observed. 

When the initial condition of the system has a projec- 
tion with one of the patterns which is much larger than 
all other projections, this overlap remains the maximal 
one and furthermore increases towards an overlap close 
to unity. Hence, the chaotic pseudo-inverse model (|4]) 
functions as an associative memory with a macroscopic 
basin of attraction (Fig. [2]). Quantitatively, the basin of 
attraction is defined as the minimal initial overlap with 
one of the patterns, and random with the other iV — 1 
directions, where the dynamics converges asymptotically 
to an average overlap greater than 0.9 with this pattern. 
Figure [5] indicates that the size of the basin of attrac- 
tion decreases with a and vanishes as a — >■ 1. The size 
of the basin of attraction also increases with the size of 
the network, N , as exemplified in Fig. [21 Similar trends 
of the basin of attractions have been observed for the 
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Hamiltonian version of the model ([T]). 

Sinee the dynamics of two different initial conditions 
inside the basin of attraction converge to an overlap close 
to unity with with the same pattern, it is important to 
verify whether the dynamics is still chaotic. The inset 
of Fig. [2] presents the Lyapunov exponent, A, measured 
from the divergence of two close initial conditions within 
an attractor for a given set of (a, e) with N — 400 and 
r = 50. Results indicate a positive Lyapunov exponent 
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(9) 



for a large range of a(l — e) < 1, as was confirmed ana- 
lytically at for a s ing le chaotic unit with infinite delayed 
feedback, r — J- oo [1^ . 

We exemplify a chaotic network with the coexistence 
of at least 2P chaotic attractors related to the unit eigen- 
vectors, i.e. the P patterns constructing Wij. Note that 
the patterns can take any real values instead of the binary 
patterns in the Ising Hamiltonian, ([1]). In the next sec- 
tion we examine the existence of other attractors besides 
the patterns and try to estimate their number. 
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III. SPIN GLASS STATES 

Starting from random initial conditions, which are typ- 
ically out of the basin of attraction and consist of overlaps 
^ 1/Vn with all patterns, result in a SG chaotic attrac- 
tor. The dynamical evolution of the P overlaps from a 
random initial condition up to 10"^ r steps is depicted in 
Fig. [5^, where each point is averaged over a window of r 
steps. The average is essential since the overlaps dramati- 
cally change from one step to another as a consequence of 
the discrete time nature of the chaotic map. The P over- 
laps varies substantially with time. In addition, cross- 
ings among the overlap amplitudes are observed during 
a transient time. 

In contrast, after lO^r steps, the average overlaps set- 
tled into a well defined order as depicted in Fig. 
A similar fixed order in the overlap amplitudes was ob- 
served up to lO^T steps. The overlaps are still not macro- 
scopic and are typically bounded (in their absolute value) 
by 3/\/N as was confirmed in simulations of network up 
to iV = 1000 units. Starting from two close initial condi- 
tions within a SG chaotic attractor, a similar Lyapunov 
exponent A was measured as for Mattis chaotic attrac- 
tors, eq. ([9|), independent of (a, e). The dynamics are 
still chaotic, but the trajectory stays (on average) in the 
neighborhood of a unique direction characterized by a 
given set of projections with the patterns. Fig. |31d. 



A. Number of Attractors 

Quantitative estimation of the freezing of the P over- 
laps requires a sophisticated method. The overlap am- 
plitudes are ranked in a decreasing order. The pattern 



FIG. 3: (color online). Spin-glass chaotic attractors are ex- 
emplified in simulations with e = 0.9 and A'' = 100. (a) The 
dynamical evolution of the overlaps with the P patterns as a 
function of time, where each point is averaged over a window 
of r steps. Crossings among the P overlaps are observed dur- 
ing the initial transient time, (b) The overlaps of panel a after 
the network settle into one of the SG chaotic attractors and 
the overlaps are constant, (c) Probability density of the Man- 
hattan distance along a trajectory while starting from a ran- 
dom initial condition. A"'". The distribution is constructed 
from the calculation of the averaged A"'" over segments of 
lO^r steps and time offset of lO'^r between consecutive seg- 
ments, (d) Probability density of the Manhattan distance 
A"'^ between the trajectories obtained from two random ini- 
tial conditions and with an offset of 10"*r. The distribution 
was constructed similar to panel c. 



with the largest overlap is labeled as P, the pattern with 
the second largest overlap as P — 1 ... and the pattern 
with the smallest overlap as 1. This ranking results in a 
vector R which is a permutation of { 1 . . -P}- Next we 
introduce the Manhattan distance llSlI 



1 P 



(10) 



between two rankings, a and /?. Note that A"'^ £ 
[0,P/2]. For a random order of with respect to Rq, 
the approximate mean value is 



1 



(11) 



where correlations among P^ — for different /i's are 
neglected. 
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Results of distribution of the Manhattan distance 
along different segments of the same trajectory a, A"'", 
when starting from a random initial condition and with 
an offset time of W^t is depicted in Fig. [5J;. The distri- 
bution is constructed from the Manhattan distance be- 
tween each pair of 100 measurements of Rq with a time 
offset between two consecutive measurements along the 
trajectory being lO'^r. Results indicate an almost com- 
plete freezing, where on the average less than a single 
change in the Manhattan order occurs after an offset of 
lO'^T steps. In order to exclude an accumulated drift in 
the Manhattan order on much larger time scales, we re- 
peated the distribution of Fig. [St with an offset of lO^'r 
and lO^r steps among successive measurements and ob- 
tained similar distributions. Figure [3Ji depicts the distri- 
bution of A"'^ between the trajectories of two random 
initial conditions and with an offset of lO^r. The most 
probable distance, A"'^, is indeed close to the approxi- 
mated one, « P/3 = 50/3 « 16, where the lower tail of 
the distribution ends approximately at 10. 

A comparison between Fig. [3J:-d for various N indi- 
cates that there is a substantial gap between the tails of 
the distributions of A"'" and A"'^ as long as TV > 30. 
This gap is the fundamental observation which enables 
us to estimate the number of different SG chaotic at- 
tractors. Assume all SG attractors have similar basin of 
attractions, one can estimate the number of attractors 
for a given N, a, e and a using the following procedure. 
We start sequentially from A'o random initial conditions 
for the network, eq. and after a time offset of at 

least lO^r the stationary ranking R is identified. In the 
event the Manhattan distances of this new identified vec- 
tor with all previously stored ranking vectors are greater 
than 1, a new SG attractor is identified. Next, its rank- 
ing vector R is added to the list of stored SG attractors. 
In order to well sample the number of SG chaotic attrac- 
tors, Nq is selected such that 7iiVo ends in previously 
identified chaotic SG attractors where 71 ss 0.1. Wc in- 
deed observed that the revisited SG attractors in this 
procedure occur almost exclusively once, which strength- 
ens our self-consistent assumption that all SG attractors 
have similar basin of attractions. At the end of this step 
of the procedure (1 — 7i)A^o distinct SG attractors are 
revealed. In the second step, we select A'^i new random 
initial conditions and calculate the fraction 72 of revis- 
ited (1— 7i)iVo recorded attractors similar to the previous 
step. The number of attractors is then estimated to be 



10 



(l-7i)iVo 



72 



(12) 



We expect this procedure to yield a lower bound for the 
number of SG attractors, since attractors with slightly 
larger basin of attractions are expected to be revealed 
with higher probability in the first step of the algorithm. 
Figure 3 depicts the estimated number of SG chaotic at- 
tractors as a function of following the above procedure 
for e = 0.5 and 0.9. The estimation was carried out in 
the range of = [30,60], since for A^ < 30 the distri- 
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FIG. 4: (color online). The estimation of the number of SG 
chaotic attractors as a function of A'' for e = 0.5 (red) and e = 
0.9 (blue). For each and e the average number of SG chaotic 
attractors and the standard deviation were derived from at 
least 10 pseudo-inverse systems. Results for the number of 
chaotic attractors for the Mackey-Glass and a = 0.5 (black) 



butions of A"^" and A"^*^ (Fig- Et-d), are significantly 
overlapped and for A^ > 60 the number of attractors is 
too large to be reliably estimated in a reasonable com- 
putational time. Results indicate that the number of SG 
chaotic attractors scales exponentially with the size of 



the system. 



with A ~ 0.19. It seems that the pref- 



actor A slightly increase with e as the local dynamics 
term in (|4]) is weakened. As a conclusion, a SG phase 
with exponentially many chaotic attractors is found. 



B. Spin-Glass Phase 

A conventional method to identify a SG phase is a 
freezing of each degree of freedom , Xi, which cannot 
be directly deduced from the freezing estimated by the 
Manhattan distance, (fTO| . To measure the freezing of Xi 
we introduce a quantity 



(A, 



(Xf) 



(13) 



where ( ) stands for an average over time steps. Eq. 
is analogous to the Edwards- Anderson order parameter, 
Qea ~ TqY^i4>i Sill- Figure [S^ indicates that along 
a trajectory inside a basin of attraction of a given pat- 
tern (blue curve) all units are almost completely frozen, 
(pi ^ 0.95. Note that for a uniform distribution of Xi 
with a given sign, e.g. Xi S [0, 1], one can verify that 
4)i = 0.75; hence (pi ^ 0.95 indicates a freezing beyond 
the preservation of the sign of each Xi . Indeed, a closer 
look at the distribution of each Xi indicates that the dis- 
tribution has considerably narrowed, e.g. mostly within 
the limited range Xi S [0.5, 1]. A typical freezing within 
a SG chaotic attractor is depicted by the red curve of 
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Fig. [5^. Almost all the units are highly frozen; how- 
ever, a slight degradation in the freezing is observed in 
comparison to the blue curve. 

The relative freezing direction between two chaotic SG 
trajectories is defined by 



(X">j| il(x^) 



(14) 



similarly to the definition of the overlap between two 
thermodynamic states of a SG system [J and q"'^ € 
[—1,1] . It is expected that the absolute overlap between 
two different SG chaotic trajectories, a ^ (3 is 

close to zero whereas the overlap between two vastly sep- 
arated time windows along a trajectory of a given SG 
chaotic trajectory, is close to unity. Figure [SJd 

such a behavior where a substantial gap between the dis- 
tributions of \q"'^\ and jg"'"] is observed, which is in 
agreement also with Fig. |3t-d. 

A crucial question is whether the SG chaotic attractors 
are closely related to the metastable states of the pseudo- 
inverse Ising spin system, Eqs. (|ll2p . with or without 
the diagonal terms. This analogy is stimulating since the 
number of metastable states for the infinite-range SG sys- 
tem also scales exponentially with TV and with A ~ 0.2 
d, 01 • We examined this question by clipping the sign of 
the frozen units in a SG attractor, 5*° = sign {{Xi)), and 
then evolved the pseudo- inverse system, eq. (P), using 
zero temperature Monte-Carlo dynamics to the nearest 
metastable state with the Hamiltonian scenario, Wu = 0, 
or with Wii 7^ and with long dynamics and obtained a 
configuration {S**}. In both cases the overlap jj J^i S^Sf 
was found to be close to zero, indicating no simple inter- 
play between metastable states of the energy surface, eq. 
([1]). and the chaotic SG attractors. The lack of correla- 
tion was also obtained in the reverse scenario where the 
chaotic dynamics, {X^}, was initiated by a metastable 
state of eq. ([T]). 



IV. CONTINUOUS TIME SYSTEMS 
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FIG. 5: (color online), (a) Freezing of each one of the N = 
100 units in a trajectory with an overlap close to unity with 
a patterns (blue-upper curve) and in a SG chaotic attractor 
(red- lower curve), (b) The overlap g"''', identifying a relative 
direction of freezing between two chaotic SG attractors (dark 
blue), and between two vastly separated (10*r steps) time 
windows along a trajectory of a given SG chaotic trajectory 
qa.a rpj^g probability density Pr{q°''^) was derived from 100 
chaotic SG attractors, whereas Pr{q"'°') was constructed from 
100 windows of length lOOr (light orange). 



and 



g{x) 



l + Ki-r)r 



(17) 



As before, N stands for the number of units in the net- 
work and the matrix Wij is given by eq. ([2]). The pa- 
rameters used in our simulations are /3 = 2,7 = l,r = 2, 
c = 10. 

Using same arguments as in ([5]), also in this case 
decreases exponentially to zero. Hence, the stability of 
the Mattis states as well as the size of the basin of at- 
tractions can be estimated as in section |TT1 Similarly, the 
number of SG chaotic attractors can be estimated using 
the same procedure used in section [III Al Results for the 
number of chaotic SG attractors for a = 0.5 are depicted 
in Fig. 21 indicating again an exponential scaling with 
the network size, TV. 



Finally we report on the extension of the exponentially 
many chaotic attractors to a network of coupled differ- 
ential equations, the Mackey-Glass (MG) equation [l3 |. 
originated from physiological control systems. This equa- 
tion was also found to have many applications in hema- 
tology, cardiology, neurology and psychiatry [lE[lB|- The 
counterpart network version of the MG equation is given 
for unit i, for instance, by 

where in the notations of eq. ^ 

F{x) = -ix{t) (16) 



V. CONCLUSIONS 

We investigated discrete and continuous time networks 
of chaotic units with delayed interactions. It has been 
shown that the networks function as an associative mem- 
ory with a macroscopic basin of attraction (Fig. [2]). 

Moreover, a SG phase which is characterized by an 
exponential number of attractors (Fig. 2]) is identified, 
where two nearby trajectories within a given attractor 
diverge from each other with a positive Lyapunov expo- 
nent ©. 

It is worthwhile to note that the slope of the expo- 
nential number of chaotic attractors is very similar for 
both the Bernoulli map and the Mackcy-Glass system, 
w 0.19 and sa 0.21, respectively, which might indicate 
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a universal behavior. Both the Bernouhi map and the is whether muhiplc attractors and a similar SG behavior 
Mackey- Glass equation are sign preserving, even in the can be found in a chaotic dynamics with the lack of sign 
generalized forms eqs. ([5]) and (|15|) . An open question preservation. 
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